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ABSTRACT 

In a previous paper, we have presented a global view of the stability of Neptune Trojan 
(NT hereafter) on inclined orbit. As the continuation of the investigation, we discuss 
in this paper the dependence of stability of NT orbits on the eccentricity. For this 
task, high-resolution dynamical maps are constructed using the results of extensive 
numerical integrations of orbits initialized on the fine grids of initial semimajor axis 
(ao) versus eccentricity (eo). The extensions of regions of stable orbits on the (ao, eo)- 
plane at different inclinations are shown. The maximum eccentricities of stable orbits 
in three most stable regions at low (0°,12°), medium (22°, 36°) and high (51°, 59°) 
inclination, are found to be 0.10, 0.12 and 0.04, respectively. The fine structures in the 
dynamical maps are described. Via the frequency analysis method, the mechanisms 
that portray the dynamical maps are revealed. The secondary resonances, concerning 
the frequency of the librating resonant angle A — As and the frequency of the quasi 2:1 
mean motion resonance (MMR for short hereafter) between Neptune and Uranus, are 
found deeply involved in the motion of NTs. Secular resonances are detected and they 
also contribute significantly to the triggering of chaos in the motion. Particularly, the 
effects of the secular resonance fg, z^ig are clarified. 

We also investigate the orbital stabilities of six observed NTs by checking the orbits 
of hundreds clones of them generated within the observing error bars. We conclude 
that four of them, except 2001 QR322 and 2005 T074, are deeply inside the stable 
region. The 2001 QR322 is in the close vicinity of the most significant secondary 
resonance. The 2005 T074 locates close to the boundary separating stable orbits from 
unstable ones, and it may be influenced by a secular resonance. 
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1 INTRODUCTION 

Besides the trans-Neptunian region (also known as "Kuiper 
belt" ) outside the orbit of Neptune and the "main belt" in 
between the orbits of Mars and Jupiter, the trojan-cloud 
around some planets is another main reservoir of minor ob- 
jects in the Solar system. Thousands of Trojan asteroids 
orbiting around the equilateral Lagrange equilibrium points 
of Jupiter have been observed and catalogued since the dis- 
covery of (588) Achilles in 1906 (Nicholson 1961). In recent 

years, six NTs were discovered! and many more of them are 
expected in future (Sheppard & Trujillo 2006). As for other 
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planets, except for Mars, of which 4 Trojans are known, 
investigations show the evidences of large unstable areas 
around the triangular Lagrange points (Marzari, Tricarico 
& Scholl 2003a; Nesvorny & Dones 2002; Dvorak, Bazso & 
Zhou 2010). 

Since the observed NTs show a wide range of orbital 
elements (especially the inclination, see Table 1) and the 
potential Trojan swarm could be much more thick than the 
one of Jupiter (Sheppard & Trujillo 2006), many studies 
have focused on this asteroid population, with interests both 
in the orbital dynamics and capturing history e.g., (Brasser 
et al. 2004; Zhou, Dvorak & Sun 2009a; Zhou, Dvorak & 
Sun 2009b; Li, Zhou & Sun 2007; Nesvorny & Vokrouhlicky 
2009; Lykawka et al 2009). 

In a previous paper (Zhou, Dvorak & Sun 2009, here- 
after Paper I), we investigated the orbital stability of 
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massless artificial asteroids around the triangular Lagrange 
points of Neptune (La and L5) in a model consisting of the 
Sun and four outer planets (Jupiter, Saturn, Uranus and 
Neptune). We constructed the "dynamical maps" on the 
plane of initial semimajor axis versus inclination (oo,io) to 
study the dependence of stability on inclination. We con- 
firmed two stable regions with inclination in the ranges 
of (0°,12°) and (22°, 36°) (Nesvorny & Dones 2002), and 
found a new stable region for highly-inclined orbits with 
io e (51°, 59°). 

To complete the investigation of orbital stability in the 
whole orbital parameter space, we present in this paper a 
detailed analysis on the dependence of orbital stability on 
the initial eccentricity. We will explain the dynamical model, 
numerical algorithm and analytical method used in this pa- 
per in Section 2. The main results will be given in Section 3, 
where the dynamical maps on the plane of (ao, eo) are plot- 
ted and described. Then in Section 4, through a frequency 
analysis method, we will discuss the mechanisms that sculp- 
ture the dynamical maps. In Section 5, the orbits of the 6 
observed NTs will be checked and compared with our pre- 
vious results. Finally, a summation and a discussion will be 
given in Section 6. 



2 MODEL AND METHOD 

Basically, the dynamical model, the numerical method and 
the analytical tools used in this paper are the same or very 
similar to the ones in Paper I. We give a brief introduction 
below. 



2.1 Dynamical model and analyzing methods 

All the computations and analysis in this paper are under 
the "outer Solar system model" in which the massless ficti- 
tious asteroids orbit around Neptune's triangular Lagrange 
points in the gravitation field of the Sun and the four outer 
giant planets. Since L4 and L5 are dynamically symmetrical 
to each other (Nesvorny & Dones 2002; Marzari, Tricarico 
& Scholl 2003a; Zhou, Dvorak &: Sun 2009a), we may inves- 
tigate only one of them without losing the generality. Below 
we will discuss only the trailing point L5. 

To simulate the orbital evolution, we use the Lie-series 
integrator (Hanslmeier & Dvorak 1984). An on-line low- 
pass filter is embedded in to remove the short-term oscil- 
lations in the outputs (Michtchenko & Ferraz-Mello 1995; 
Michtchenko et al. 2002). The spectra of each orbit are then 
calculated from these filtered outputs using the Fast Fourier 
Transform (FFT) method. The number of peaks above a 
given level in a spectrum is counted. This so-called "Spec- 
tral Number" (SN) can be regarded as an indicator of the 
regularity of an orbit (Michtchenko & Ferraz-Mello 1995; 
Michtchenko et al. 2002; Ferraz-Mello et al. 2005), and it is 
used to construct the dynamical maps on the initial orbital 
parameter plane. The frequencies obtained from the FFT 
procedure are also carefully analyzed to determine the mech- 
anisms that portrait the features in the dynamical maps. For 
details of the above mentioned methods and technologies, 
please refer to Paper I. 
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Figure 1. The mean value (solid curve) and the libration ampli- 
tude (error bars) of a versus initial <tq. These 101 orbits have 
the same orbital elements except their mean anomalies differ 
from each other so that they have different initial resonant ar- 
gument ao. The initial eccentricity and inclination of these or- 
bits are eo = 0.1, io = 0°. The minimum of libration amplitude 
(Act = 6.97°) is at a c = -64°. 



2.2 Initial conditions and the libration center 

In order to get representative initial conditions for numeri- 
cal simulations, we follow a similar scheme as Nesvorny and 
Dones (2002) did. The initial orbital elements of fictitious 
asteroids are taken from a grid on the plane (ao, eo), where 
101 values of the semimajor axes ao are chosen uniformly 
between 29.9 and 30.5 AU, and the eccentricities eo vary 
up from 0.0 to 0.5 with an increment of 0.005. For all or- 
bits on the grid, the initial inclinations io are fixed at a 
given value. The angular elements of the orbits, including 
the ascending node Qo, the perihelion argument ojo and the 
mean anomaly Mo, are arranged in such a way that the res- 
onant angle a (a — A — Ag, the difference between the mean 
longitudes of Neptune As and the fictitious trojan A) is at 
the center of tadpole orbit (<r c ~ —60°). Specifically, we set 
fto = fis, 0J0 = u» - 60° and M = M 8 + a c + 60°. 

It is known that the libration center a c changes with 
the eccentricity and inclination. In the restricted three body 
model consisting of the Sun, Neptune and the Trojan, the 
value of a c at a given eccentricity and inclination can be 
determined by locating the extremes of the averaged Hamil- 
tonian (energy) on the (a, a) -plane (Nesvorny et al. 2002). 
In the outer solar system model adopted in this paper, be- 
cause of the perturbations from other planets, there could 
be no such a center with zero libration amplitude of a. How- 
ever, we can still locate the position at which a librates with 
the smallest amplitude. A series of orbits are integrated and 
the one with minimum libration amplitude is "defined" as 
the libration center. 

In this series of orbits, the initial semimajor axes are 
chosen as the same value as Neptune's, ao = as and eo,io 
are given the specific values. The initial angular elements 
are given the representative value as mentioned above, wo = 
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Figure 2. The variation of libration center (<r c ) with respect to 
eccentricity. From bottom to top, the curve represents the case 
for initial inclination from 0° to 60° , respectively. 



oj 8 - 60°, Qo = fis, Mo = Ms + 60° + a. We test a in the 
range ( — 100°, —50°). The orbits are then integrated up to 
~ 3.3 x 10 5 yr, which is about 4 times longer than the typical 
long-term variation period of a (Murray & Dermott 1999) . 
We determine the "practical" libration center from the be- 
haviour of the mean value and the librating amplitude of a. 
An example with initial eccentricity eo = 0.1 and inclination 
id — 0° is shown in Fig. 1. 

It is interesting to notice in Fig. 1 that around the libra- 
tion center the libration amplitude increases with the differ- 
ence of the initial ao from the center a c and both of the 
mean value and oscillation amplitude change continuously, 
reflecting the regularity of motion around the libration cen- 
ter. When the initial <7o is far away from the center, e.g. 
(To < —90°, the libration amplitude becomes large and the 
mean value changes randomly, indicating the chaotic char- 
acter of motion. 

Varying eo,io and repeating the above calculations, we 
find the variation of libration center with respect to eo and 
io. A summation is plotted in Fig. 2 for orbits with initial 
conditions of eo = — 0.30 and io = 0° — 60°. We ignore 
higher eccentricity and inclination because the highly eccen- 
tric and/or inclined orbits are unstable and therefore are of 
less interest. In fact, at very high eccentricity, the libration 
center may "bifurcate" to two centers. But this phenomenon 
is out of the range of current paper. 

As shown in Fig. 2, the libration center of the tadpole 
orbit shifts farther away from Neptune as the eccentricity 
of the orbit increases. For the coplanar orbit (io = 0°), the 
libration center is exactly at the equilateral triangular point 
with a c = —60° when the eccentricity is 0, and it changes 
continuously to a c ~ —71° when the eccentricity increases 
to 0.25. The o c also changes a little with respect to the in- 
clination, and compared with coplanar orbits, the variation 
of <r c with eccentricity for highly inclined orbits becomes 
smaller. 

Only after knowing the libration center cr c (eo, io) at dif- 
ferent eccentricity (eo) and inclination (jo), we can select the 



proper initial conditions to construct the representative dy- 
namical maps. Namely, for a given inclination io, a dynam- 
ical map on the (ao, eo)- plane is constructed by analyzing 
the stability of each orbit on this plane and with other initial 
angular elements given by fio = ^8,^0 = <^8 — 60°, Mo = 
M 8 +60° + <7 c (e ,io). 



3 DYNAMICAL MAPS 

We have investigated in Paper I the dependence of orbital 
stability of NT on the inclination. In this paper, we focus 
on the dependence of stability on the eccentricity. For this 
sake, the initial inclinations io is fixed at given values and 
we study the orbits of the fictitious trojans initialized on 
the grid on the (ao, eo)- plane with other orbital elements 
determined in a way described in last section. Using the 
"spectra number" (SN) as the indicator of the regularity of 
an orbit as in Paper I, we construct dynamical maps on the 
(a ,e )- plane. 

We present in Fig. 3 the dynamical maps for io = 
0°, 10°, 20°, 30°, 40° and i = 50°. The color in these dy- 
namical maps indicates the orbital regularity. The green 
represents the most regular orbit, the red is on the edge 
of chaotic motion and the white indicates that the orbit 
does not survive on the trojan-like orbit in our numerical 
simulation (which lasts about 34Myr). According to our ex- 
periences in Paper I, orbits with SN> 50 (color blue) are not 
to be expected to survive in the Solar system age. Surely it 
is impossible and also unnecessary to calculate each slice at 
every possible initial inclination. These slices in Fig. 3 are 
enough to show the variation of dynamical map with re- 
spect to the initial inclination. We also discard those cases 
in which most of the trojan orbits are unstable and there- 
fore less interested. For instance, due to the depletion effects 
arising from the z/8 secular resonance and Kozai resonance, 
the motion of trojans with io ~ 45° and io > 60° are un- 
stable (see Paper I), so that the dynamical maps are plain 
and not shown here. However, for the sake of providing more 
information about the trojan motion at other inclinations, 
other three dynamical maps for io = 5°, 35° and 55° are 
used in Figs. 5, 6 & 7 as examples to present our frequency 
analysis in next section. 

The first conclusion drawn from these dynamical 
maps may be that the most stable motion (color green) 
does not happen at any inclination, but only at io = 
0°, 5°, 10°, 30°, 35° and 55°. This is consistent with the find- 
ing in Paper I of three most stable regions in inclination: 
Region A (i = 0° - 12°), B (i = 22° - 36°) and C 
(io = 51° -59°). 

A second conclusion is that the initial eccentricity eo 
for stable orbits must be small. Although for some initial 
inclinations (e.g. io = 30°, 35°) some orbits with eo as large 
as 0.25 survive in our integrations and may survive even for 
much longer time, their SN's (color red) indicate that most 
probably they will not remain on the trojan-like orbits in 
the Solar system age. These orbits with high eo may suf- 
fer perturbations through resonances whose effects appear 
significantly only in a timespan much longer than our inte- 
gration time. For the most stable orbits, the largest initial 
eccentricities are around 0.08 in Region A, 0.07 in Region B 
and 0.02 in Region C. However, if we use a lenient criterion, 
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Figure 3. Dynamical maps on (ao, eo)— plane, for initial inclinations 0°, 10°, 20°, 30°, 40° and 50° (left to right, top to bottom). The 
abscissae and ordinates are initial semimajor axes and eccentricities. Note the ordinates may have different scales. The contour curves 
represent the libration amplitude of resonant angle ct = A — Ag. From inside out are contours for Act = 10° to Act = 70° with an 
increment of 10°. Particularly, the red curves are for Act = 30° and the thick ones represent Act = 60°. 
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namely, regarding all the orbits with SN< 50 as survival 
ones (in the Solar system age), the permitted eccentricity 
extends to 0.10, 0.12 and 0.04 in Region A, B and C respec- 
tively. Out of these three most stable regions in inclination, 
there are also some narrow areas in which the motion is sta- 
ble. For example, on the slice of io = 20° in Fig. 3, we see a 
few stable orbits with eo < 0.05. 

The contour of libration amplitude of the resonant an- 
gle Act = CT max — cr m i n is also calculated and over-plotted 
on the dynamical maps. Clearly the relation between Act 
and the deviation of semimajor axis from its mean value 
a -a c (Erdi 1988; Milani 1993; Zhou, Dvorak & Sun 2009a) 
is obeyed. The Act increases monotonically with the distance 
of initial semimajor axis ao from its value at libration center 
(a c ~ 30.22 AU). The maximum Act is around Act ~ 70° for 
stable orbits with low (Region A) and medium inclination 
(Region B), but shrinks a little for high inclination (Region 
C) to A ~ 50°. At given io and ao, the libration ampli- 
tude increases only very little as the initial eccentricity in- 
creases. At high eccentricity, the Act may still be small, but 
the large SN implies a chaotic motion. Therefore, the small 
libration amplitude manifested in our integration time may 
increase in a distant future. We also note that there is a 
lack of small Act at some inclination values. For example, 
at i = 0°, 10° and 50° all trojan orbits have Act > 10°. 
Last but not least, thanks to the calculation of libration 
center ct c and proper selection of initial conditions, both the 
dynamical maps and the libration amplitude contours have 
the bilaterally symmetrical features. Otherwise, if we just 
simply set ct c = —60°, these figures would be distorted. 

Contrary to the profiles of the libration amplitude con- 
tours that are more or less plain, the dynamical maps bear 
plenty of fine structures. Most distinctly, there are verti- 
cal gaps of irregular motion on slices of small inclination 
(io = 0° , 5° , 10° ) , and a "spike" of relatively regular orbits 
with small Act can be clearly seen on slices of large incli- 
nation (io = 40°, 50°, 55°). On one hand, these fine struc- 
tures reveal the peculiar behaviour of orbits in these regions; 
on the other hand, they supply hints of the mechanisms 
taking effects and portraying the dynamical maps in the 
corresponding parameter space. We will try in next section 
to locate different kinds of resonances on the (ao,eo)- plane 
and compare them with the fine structures in the dynamical 
maps. 

We would like to notify that the dynamical maps are 
plotted based on the osculating orbital elements in our pa- 
pers. As we have shown in Paper I, the dynamical map would 
change if a different epoch was chosen. In a different epoch, 
the features of the dynamical map are very well kept, and 
the change is nothing more than a shift in the semimajor 
axis. If the proper libration amplitude and proper eccentric- 
ity, instead of the osculating semimajor axis and osculating 
eccentricity, are used to define the representative plane for 
the dynamical map, we may obtain the epoch-free dynami- 
cal maps as in (Marzari, Tricarico & Scholl 2003a; Marzari, 
Tricarico & Scholl 2003b; Dvorak et al. 2007). But the os- 
culating elements have the merit of being able to compare 
with the observing data directly. To do this, all we need is 
to transform the observing data to the corresponding epoch, 
as we will do in Section 5 of this paper. 

We have calculated many other dynamical maps with 
different initial inclinations, of which three maps for io = 




20 30 40 

Initial Inclination (i ) 

Figure 4. The number of regular orbits versus the initial in- 
clination. From bottom up, the solid squares, solid circles, open 
squares, open circles and stars indicate the numbers of orbits 
whose SN's are smaller than 25, 45, 65, 85 and 105. Remember 
all those orbits with original SN > 100 and 29.9 au < a < 30.5 au 
are assigned a new SN of 100 (please refer to Paper I for details). 



5°, 35° and io = 55° are going to be shown in next section. 
Instead of presenting all of them here, we count the number 
of regular orbits on the dynamical map of each selected in- 
clination value. Since for each map we have calculated 5151 
orbits on a 101 x 51 grid of (fflo,eo), the number of orbits 
with the SN smaller than a definite value is an indicator 
of the relative regularity of orbits. The results are summa- 
rized in Fig. 4, where we show the numbers of orbits with 
SN smaller than 25, 45, 65, 85 and 105. According to our 
rule of assigning SN, all asteroids that escaped from thel:l 
MMR (judged by their mean semimajor axes during the or- 
bital integrations) have been designated an SN of 110. Thus 
the stars in Fig. 4 imply the surviving probability of artifi- 
cial Trojans out of those totally 5151 ones in our integration 
time. 

The survival number of artificial Trojans (stars in 
Fig. 4) increases with the initial inclination before it reaches 
its maxima at io = 30°. Then it drops down to a mini- 
mum at io = 45° due to the depletion by the u& secular 
resonance, and another maximum is reached at io = 50°. It 
ends at higher inclination where the Kozai resonance domi- 
nates (Paper I). If adopting the survival probability from our 
integrations as the stability indicator, we would conjecture 
that more NTs could be found on high-inclined orbits with 
inclination around 30° and 50°. But if we focus on orbits 
with SN < 25 and/or SN < 45, we would argue that NTs 
may mainly concentrate in a low-inclination region around 
io ~ 5° and a medium-inclination region around io ~ 30°. 
This has been partly proven to be true by the observations, 
namely, among the six observed NTs, 4 are found in the 
former region and the rest 2 in the latter. In addition, judg- 
ing from Fig. 4, a high-inclination region around io ~ 55° 
also could host some NTs. However, we would note that the 
stability of the high-inclined orbits may be determined by 
the dynamical analysis under current configuration of the 



© 2009 RAS, MNRAS 000, 1- 



6 Zhou, Dvorak & Sun 



Solar system, but the existence of real NTs on the specific 
orbits depends also on the evolution of planetary orbits in 
the early stage of the Solar system (Kortenkamp, Malho- 
tra & Michtchenko 2004; Li, Zhou & Sun 2007; Nesvorny & 
Vokrouhlicky 2009). 



4 FREQUENCY ANALYSIS 

4.1 Basic ideas 

As we have mentioned above, some secondary resonances 
and secular resonances are responsible for the forming of the 
fine structures in the dynamical maps. To understand the 
dynamics of NTs clearly, we need to find the resonances as- 
sociated with these structures. One natural and simple way 
to determine a resonance is to check the behaviour of the 
corresponding critical resonant angle, whose librating char- 
acter is an evident sign of being inside a resonance. But, as 
mentioned in Robutel & Gabern (2006), this method is valid 
only for low order MMRs or secular resonances of second or- 
der. It doesn't work when the resonance involves more than 
two degrees of freedom or when it is of high order. In fact, 
a typical power spectrum of a complex motion generally is 
a composition of an amount of terms with comparable am- 
plitudes and different frequencies, this means, the regular 
variation of a critical angle of a high-order resonance may 
be enshrouded by other high-amplitude terms even when the 
motion is deep inside the resonance. Of course, the proper 
critical angle can still be well defined if we can remove the 
irrelevant high-amplitude terms by (for example) reducing 
the corresponding Hamiltonian to a normal form (Robutel 
& Gabern 2006). To define such critical angles is out of the 
scope of this paper, so that we will not try to check whether 
the resonant angle is librating or circulating when we dis- 
cuss a resonance. Instead, two techniques are applied in this 
paper. One is the "dynamical spectrum" and the other is 
the "empirical formulation". These two methods have been 
used in Paper I, and we will give a very brief introduction 
below. 

For each NT orbit in our simulation, we calculated the 
power spectra of a cos a, e cos vj and icosfl These spectra 
generally are complex compositions of the forced frequen- 
cies, free frequencies, their harmonics and combinations of 
above terms. From a spectrum, the most significant peaks 
are selected and their frequencies are recorded. 

We vary the orbital parameter and calculate these lead- 
ing frequencies at each parameter value. By plotting all these 
frequencies on one figure (called "dynamical spectrum" as in 
Paper I), we may find how the dominating frequencies of the 
motion change with respect to the parameter. Knowing the 
secular frequencies of the outer solar system, i.e. the apsidal 
and nodal precession frequencies of Jupiter, Saturn, Uranus 
and Neptune (g$, <?6,<?7, <?8 and S5, s@, sj, ss), the forced fre- 
quencies can be easily recognized. Since the proper frequen- 
cies of the motion generally change continuously with the 
orbital parameter, they can also be recognized. In a dynam- 
ical spectrum, along with the varying parameter, a contin- 
uously varying frequency may meet and cross another fre- 
quency, where a resonance associated with these frequencies 
happens. 

Thus a dynamical spectrum can be used to detect and 



locate the resonances in the motion, and more importantly, 
it helps us determine the proper libration and precession fre- 
quencies of the motion. The proper frequencies of the NT's 
motion include the libration frequency of the critical angle 
of the 1:1 MMR, / CT , the frequencies of the apsidal preces- 
sion g and of the nodal precession s. They can be derived 
from the dynamical spectra of the variables a cos a, e cos vj 
and i cos Q respectively. With the data of the proper fre- 
quencies at different parameter values in hand, we can com- 
pute the empirical formulations for the proper frequencies 
with respect to orbital parameters (e.g. the semimajor axis, 
inclination and eccentricity) on the whole parameter space. 
Finally, with these empirical formulations and those already- 
known frequencies, we may calculate all kinds of equations 
of frequencies that define the resonances in the parameter 
space. 

Since the dynamical spectra are generally overloaded 
with details (see for example Figs. 9&10 in Paper I), we 
would rather not show and describe them here. Instead, we 
will give directly the empirical formulae of f a ,g,s derived 
from the dynamical spectra, and then with the formulae we 
figure out the resonances revealed by the dynamical spectra 
and manifested in the dynamical maps. 

To obtain a global view of the dynamics of NTs at differ- 
ent inclinations, we analyze three typical dynamical maps on 
the (ao, eo)— plane with io = 5°, 35° and 55°, respectively. 
These selected slices are inside the most stable regions Re- 
gion A, B and C. 

4.2 For i = 5° 

Adopting the similar quadratic formula (Milani 1994; 
Marzari, Tricarico & Scholl 2003b) as in Paper I, and setting 
8 = ao — 30.215, we derive from our calculations the empiri- 
cal expressions of fa,g,s on the (ao, eo)— plane for the slice 
io = 5°, which are: 

/ CT [10~ 4 27r/yr] = 1.1350 - 2.17835 2 - 0.046311e 2 



+ 2.0566o - 27.500e 4 + 28.700o e , 
5 [10~ 6 27r/yr] = 1.5783 + 6.9896<5 2 + 3.3556e 2 

- 72.941o 4 + 174. 70e 4 - 110.97<5 2 e 2 , 
s[10~ 7 27r/yr] = 5.1205 + 21.435<5 2 + 19.102e 2 

- 54.217o 4 + 3567.4e 4 - 416.73o 2 e 2 . 



(1) 



(2) 



(3) 



Similar calculations are performed for slices of io = 35° , 55° 
and will be presented below. Using these formulae, all the 
locations of different kinds of commensurabilities between 
frequencies can be computed. However, before doing this, 
let us recall what kinds of resonances may happen in NTs' 
motion. 

In the study of Jupiter Trojans, Robutel & Gabern 
(2006) defined four families of resonances as follow: 



I : pf a - n 5 + qg + q 5 g 5 + qegt, = 0; 

II : f a + 5/ 5; 2 +pg + P555 + P6P6 = 0; 
III: qs + q 6 s,j + p 5 g 5 + p 6 g 6 = 0; 

IV : pg + / 5 :2 + psg-> + pega = 0. 



(4) 



In above equations, 715 is the mean motion of Jupiter, /s : 2 is 
the frequency of the quasi 5:2 resonance between Jupiter 
and Saturn (the Great Inequality), and all the integers 
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p,q,P5,pe,qs and q% should be chosen in a way such that 
the d'Alembert rule is fulfilled. 

It is reasonable to postulate that similar resonances 
are involved in the case for NTs' motion. Among the 
already-known frequencies, except for those secular frequen- 
cies p5, <?6, <?7,<?8 and s 5 ,se,S7,ss, we find the frequency of 
the quasi 2:1 MMR between Neptune and Uranus, that is, 
the frequency of the angle 2As — A7 (denoted by f2.i here- 
after), plays an important role in the motion of NT. On the 
contrary, the perturbing effects from the Great Inequality 
is (relatively) ignorable, due to its far distance and perhaps 
also due to its high frequency (~ IO -3 27r/yr). 

In addition, we notice that the mean motion of Neptune 
ns(~ 6.068 x IO -3 2-7r/yr) is much higher than the libration 
frequency /„(~ 1 x 10 -4 2-7r/yr). This means, if a resonance 
similar to Family-I happened, the integer p would be as large 
as ~ 60. Moreover, / 2: i = 2.36064 x 10~ 4 27r/yr is at least 
2 orders of magnitude larger than the precession rate g, s 
(and also other precession frequencies related to planets). 
Therefore, we may ignore those resonances as Family-I and 
Family-IV in NT's motion, because their order must be quite 
high. Instead, only two types of resonances are discussed in 
this paper: 



0.12 



S : 

c 



pg + qs + Yfj =B (Pi9i + QjSj) = 0, 

hfa + kf 2:1 +pg + qs + Yf j= 5(Pi9j 



■qjSj) = 0, 



(•») 



where h,k,p,q,pj,qj are integers. The "S" stands for "Sec- 
ular", indicating they are secular resonances in the usual 
sense. While the "C" for "Combined" indicates that these 
resonances are combinations of mean motions and secular 
frequencies. The d'Alembert rule for S-type resonance is 
P+3+Z^ =B (Pj+3j) = °. and it's k+p+q+Y? j=s (pj+qj) = 
for C-type. In both equations, (q + ^~J qj) must be even 
so that the symmetry of inclination with respect to the ref- 
erence plane can be guaranteed (Murray & Dermott 1999). 

By checking carefully the details in the dynamical spec- 
tra and using the empirical expressions in Eqs. (l)-(3), we 
recognize some important resonances taking effects on the 
slice of io — 5°, whose locations are calculated and then 
plotted over the corresponding dynamical map as shown in 
Fig. 5. 

The most distinguishable structure in the dynamical 
map of io = 5° (and also other slices at low inclination, e.g. 
io = 0° , 10° in Fig. 3) is the blue vertical "stripes" indicating 
less stable motion. We find that they arise from the effects 
of some C-type resonances, among which three major ones 
are shown in Fig. 5. The dashed curves labeled by CI, C2 
and C3 are locations of the following resonances: 



CI 
C2 
C3 



2/„-/2:l+flB=0, 

4/ CT -2/ 2:1 +.g 6 + 37 = 0, 

6/cr - 3/ 2 :l + #5 +36+37 



= 0. 



(6) 



Apparently, their locations match the positions of vertical 
structures very well. We need to stress here that each res- 
onance listed above (and below too) is just a representa- 
tive of a bunch of resonances, characterized by higher or- 
der combinations of integers p,q,pj,qj in Eq. (5). For ex- 
ample, beside the resonance CI that reaches the abscissa 
at ao = 30.044 au and 30.392 au, we can easily find another 
two resonances as CI' : 2/ CT — fc-.i — g + gs + g§ = and 
CI" : 2/ CT — /2:i +g — <?8 + <?6 = 0. And in fact, the resonance 
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Figure 5. The locations of resonances plotted over the dynami- 
cal map for the case with initial inclination io = 5°. The C-type 
resonances are indicated by dashed curves while the solid curves 
stand for the S-type resonances, particularly, the secular reso- 
nance 1/18 I s plotted by the thin solid curve. See text to find the 
meaning of labels CI, C2, SI, S2 and S3. 



CI itself is favourable for the stability of Trojan orbit. As 
we will see in next Section, the asteroid 2001 QR322 is deep 
inside such a resonance and its orbit is stable over the Solar 
system age. On the contrary, the resonances CI' and CI" are 
responsible for the unstable gaps in both sides of resonance 
CI, especially the apparent ones around ao = 30.050 au and 
30.375 au approximately. 

Another mechanism that may contribute to the forma- 
tion of the vertical structures is the secular nodal resonance 
1/I8 defined by s — sg = 0. Its location is also plotted in Fig. 5, 
as well as other three S-type resonances: 



SI 
S2 
S3 



3ff - 37 - gs = 0, 

g - 2s + g s - 2s$ + 2s 5 

3s - 2g 8 - Ss - 0. 



(7) 



We notice inconsistent statements about the resonance 
1/18 in the literatures. Nesvorny and Jones (2002) found in 
the dynamical map a less stable region with libration am- 
plitude between 40° and 50° (Fig. 2(d) therein), and they 
argued that this arose from the effects of the i/is resonance. 
While Brasser et al. (2004) showed that the clones of aster- 
oid 2001 QR322 "remained deep inside" the vi& resonance 
are the most stable ones ("have the longest e- folding time"). 
Marzari and Scholl (2003a) argued that this secular reso- 
nance were "not strong enough to cause a short term insta- 
bility". Our results support the same conclusion. 

In Fig. 5, no apparent details in the dynamical map can 
be found around the location of i/is resonance. In fact, we 
found in our simulations that nearly all orbits with low in- 
clination (io < 5°) and not too large libration amplitude 
are involved in the vi$ resonance, but this resonance does 
not show any prominent effects in protecting or destroying 
their orbital stabilities. The less stable region mentioned in 
Nesvorny and Jones (2002) in fact corresponds to the blue 
stripes at ao = 30.12 au and 30.30 au in the dynamical map 
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of e'o = 0° in Fig. 3, and it is clearly related to the C2 reso- 
nance in Eq. (6), but not the i>\% resonance. 

Other secular resonances, e.g. the SI, S2 and S3 in 
Eq. (7), however may have more evident dynamical effects. 
The resonances SI and S2 define the edge of the survival 
region, and the curvature of S2 in the eccentricity range of 
0.07 — 0.09 matches the outline of the stable region. The 
border between the stable region and chaotic region makes 
a V-shape blue area from eo = 0.05 to eo = 0.09, in which 
several crossovers or close encounters between resonances 
S2, S3, C2, C3 and vis happen. Thus it seems reasonable to 
argue that the border is made of a net of these resonances. 

Another observing from Fig. 5 and other dynamical 
maps at different initial inclinations in Fig. 3 is that the 
Cl-type resonance with h = 2, k = — 1 in Eq. (5) has quite 
promising dynamical effects in all initial inclination values, 
but the influences from C2 and C3 decreases as io increases. 
The blue gaps corresponding to C2 and C3 become more 
and more vague as io gets bigger. We will see this trend 
clearly in next example of io = 35° . 

4.3 For i = 35° 

Similarly, setting 5 — ao — 30.220, we obtain the best fits of 
the proper frequencies formulae for the slice io = 35°: 

/< T [10" 4 2 7 r/yr] = 1.0109 - 1.59235 2 + 0.45999e 2 

- 2.18605 4 - 29.100e 4 - 5.13705V, 
s[10 -7 27r/yr] = 8.3905 - 4.92575 2 - 6.113e 2 

- 241.095 4 - 75.306e 4 - 347.565 2 e 2 , 
s[10~ 7 27r/yr] = 2.8989 + 15.9205 2 + 7.1909e 2 



62.5025 4 + 161. 75e 4 + 368.27o 2 e 2 . 



(8) 

(9) 

(10) 



CI 


2/cr- 


- f 2:1 - g + 2g 6 = 0, 


C2 


2/ CT - 


- fc-.i - s + 2s 6 - s 8 + ga = 


C3 


2/„- 


-/ 2 :l+ ff + 2 36 -2(7 5 =0, 


C4 


2/„- 


- / 2: i +2g + 2 9(i - 3 55 = 0, 


C5 


4/ CT - 


- 2/2:i +3(?6 +gs-2g 8 = 0. 



Then by carefully studying the dynamical spectra, we detect 
some resonances involved in the motion. Using the above 
formulae, we calculate their locations and plot them over 
the corresponding dynamical map in Fig. 6. 

As in the previous case, we find in our calculations some 
C-type resonances involved in the motion of orbits with 
io = 35°, of which several typical ones are listed below and 
plotted in Fig. 6. 



(11) 



The low-order C-type resonance like the CI defined in 
Eq. (6) does not appear in the region around io = 35°. 
Among the resonances in Eq. (11), only the one with the 
lowest order (CI) has the distinguishable dynamical effects, 
shaping the sharp edge of stable region at low eccentricity 
(eo < 0.05). The C-type resonance with h — 4, k = —2 in 
Eq. (5), here C5 in Eq. (11) contributes to the formation of 
the vertical structure at the center (around ao = 30.22 au). 
Another interesting resonance is C2, in which the nodal 
precessions (s, sg and se) are involved. In fact, as the in- 
clination of Trojan orbit increases, it is natural to expect 
that the nodal precession is becoming more important. We 
would like to note that there are less regular indentations 
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Figure 6. The same as Fig. 5 but for io = 35°. Except the reso- 
nances labeled by C1-C5 (dashed curves) and S1-S4 (solid curves. 
See text for their meanings), the locations of us,uis resonances 
(solid curves in red) and Kozai resonance (thin solid curve in 
black) have also been plotted. 



corresponding to this resonance in the dynamical map at 
ao = 30.08, 30.36 au and eo > 0.015. On the other hand, al- 
though the resonances C3, C4 in Eq. (11) can be detected in 
the dynamical spectra, no apparent details in the dynamical 
map can be found related to these resonances. 

The influences from the C-type resonances are less im- 
portant in slice io = 35° than in slice io = 5°. Relatively, 
the effects from the S-type resonances are remarkable. Be- 
side the v$ , i/~L$ and Kozai resonances that locate outside the 
survival region as show in Fig. 6, several other S-type reso- 
nances as below are determined and illustrated. 



SI 
S2 

S3 

S4 



g + s - gs - sa — 0, 

g - 2s + g s - s 8 + s 5 = 0, 
g + 2s - g 6 + s 6 - 3s 5 = 0, 
g - 2s + g a + sg + gj - g 5 - s 5 



(12) 



0. 



Clearly the S2 resonance combined with S3 defines the 
boundary of the survival region, and the resonance S4 estab- 
lishes the edge of the stable region up to eo = 0.12. Again, 
on the top part of the dynamical map, the C-type and S- 
type resonances gather, resulting in overlaps among them 
and leading to chaotic motion. 

4.4 For i = 55° 

Setting S = ao — 30.223, the formulae of the proper frequen- 
cies for the slice io = 55° are: 

/ CT [10" 4 2yr/yr] = 0.88663 - 1.19985 2 - 2.0233e 2 

- 7.4512o 4 + 353.60e 4 - 
g[10~ 7 2yr/yr] = 2.7724 - 46.683<5 2 

- 1588.8o 4 + 666.87e 4 - 
s[10~ 7 27r/yr] = 1.5643 - 46.165o 2 - 52.989e 2 

+ 2168.45 4 + 5116. 5e 4 + 18800<5 2 e 2 . (15) 



106.005 e, 


(13) 


-71.157e 2 




3776.5(5 2 e 2 , 


(14) 
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Figure 7. The same as Fig. 5 but for io = 55° . Except the res- 
onances labeled by CI, C2 (dashed) and S1-S3 (solid. See text 
for the meanings), the locations of vig resonances (red lines) and 
Kozai resonance (thin solid curve in black) have also been plotted. 



The involved important resonances are determined and 
illustrated in Fig. 7, where the following C-type and S-type 
resonances are plotted: 



CI 
C2 

SI 

S2 
S3 



2/a - h-.i - 2g + 3g 6 = 0, 

2/ CT - J2:i - g — 8 + 3s 6 = 0, 

2s - 2g 5 - 2s s + ge - s a + g 7 + s 7 = 0, 

s - 2g 7 + 2s 7 - Ss — 0, 

s - g s — s 8 — g 7 + gs + s 5 = 0. 



(16) 



Although the most part of the CI lines in Fig. 7 are out of the 
survival region, we see the fine structures (narrow unstable 
gaps) around oo = 30.06, 30.36 au and eo > 0.002 are clearly 
related to this resonance. Inversely, the stability of orbits 
benefits from the C2 resonance, as we see the concentration 
of stable orbits around the C2 lines, especially in the regions 
a = 30.13, 30.31 au and e < 0.002. 

Comparatively, the dynamical effects from S-type reso- 
nances are more prominent. As shown in Fig. 7, the SI res- 
onance sharply defines the out edge of the survival region; 
the lines of S2 coincide with the boundary separating the 
chaotic orbits from the relatively regular ones; and the S3 
curves surround the most stable orbits in this slice. 

Considering the high inclination of this slice at io = 55° , 
it is not a surprise to see that the nodal precession s is 
involved more deeply in these important resonances com- 
pared to the low inclination cases (e.g., resonances listed in 
Eqs.(6), (7), (11) & (12)). 

Meanwhile, the fig resonance (red lines) and Kozai res- 
onance (thin black line) are calculated and shown in Fig. 7. 
The ^i8 is outside the survival region, that is, it does not 
play an important role in the long-term stability of orbits. If 
an orbit was trapped in Kozai resonance, its inclination i and 
eccentricity e may experience joint variations: i decreases as 
e increases, or vice versa (Kozai 1962). This variation may 
lead to a highly eccentric orbit and result in orbit escaping. 
Therefore, the boundary between the surviving and escap- 



ing orbits is determined cooperatively by the SI and Kozai 
resonances. 



4.5 A summary 

At different inclinations, different resonances participate in 
featuring the dynamical maps. Basically, the secondary reso- 
nances (C-type) related to the libration of the 1:1 MMR be- 
tween Neptune and its Trojans and to the quasi 2:1 MMR 
between Neptune and Uranus, as well as the secular reso- 
nances (S-type) related to the nodal and/or apsidal proces- 
sions, are the most important dynamical mechanisms that 
form the fine details in the dynamical maps and determine 
the complicated orbital behavior of NTs. Compared to the 
cases in low inclination, those highly inclined orbits are more 
deeply affected by the nodal-type resonances. In all cases, 
different resonances with different orders make a resonance 
net in the parameter plane of initial conditions. Thus an 
orbit may diffuse slowly along the net and migrate on the 
parameter space from one configuration to another. Of par- 
ticular interest is the slow diffusion from low inclination or- 
bits to high inclination orbits, which can be used to explain 
the origins of the observed NTs with high-inclinations, i.e. 
2005 TN53 and 2007 VL305 as listed in Table 1. 

We only reported in previous subsections some reso- 
nances detected by the dynamical spectrum method. Surely, 
many more resonances, particularly those with high orders, 
have not been revealed in our current paper, and they may 
also play important roles in characterizing the dynamical 
behaviours of Trojans' orbits, which are our studying topic 
in future. 



5 ORBITS OF OBSERVED NEPTUNE 
TROJANS 

Up to now, six NTs have been observed, and their orbital 
stabilities have been analyzed individually, e.g. in (Marzari, 
Tricarico & Scholl 2003a; Brasser et al. 2004; Li, Zhou & Sun 
2007) . Since we have presented global view of the dynamics 
of NTs, we can locate the observed NTs on the correspond- 
ing dynamical maps and draw conclusions directly from the 
comparison. 

We have adopted the configuration of the outer solar 
system at epoch of August 1, 1993 (JD2449200.5) as the 
initial conditions to simulate the orbital evolutions of four 
planets and artificial Trojans, the dynamical maps (Fig. 3 
in this paper and also Figs. 2 & 3 in Paper I) are based on 
the instantaneous orbital elements at that moment. But the 
orbits of observed NTs are given in recent epoch, e.g. the 
ones listed in Table 1 obtained from the AstDyS website^ 
are given at epoch JD2455000.5. To locate their orbits on the 
dynamical maps, we need to transfer their orbital elements 
at this moment to the epoch when the planets' orbits are 
adopted. The transferred orbits are calculated and listed in 
Table 1. 

Because all the six objects are on nearly circular or- 
bits with eccentricities smaller than 0.07, we can approxi- 
mately locate these orbits on the (ao,io)— plane in Paper I 

■f http://http://hamilton. dm. unipi.it/astdys 
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Table 1. Orbits of observed Neptune's Trojans at epochs JD=2455000.5 (June 14, 2009, six columns in the left half) and JD=2449200.5 
(August 1, 1993, six columns in the right half) in the reference frame of J2000.0. The mean anomaly M, perihelion argument u>, ascending 
node f! and inclination i are in degrees. All the angle values have been rounded to one decimal place for simplicity. The last column Act 
gives the approximations of the libration amplitude (in degrees) of a, derived from Fig. 2 in Paper I. 



Designation 



M 



n 



a(AU) 



M 







a(AU) 



2001 QR322 


57.31 


162.5 


151.6 


1.3 


0.0316 


30.336 


— > 


30.06 


155.3 


151.6 


1.3 


0.0290 


30.208 


77 


2004 UP10 


343.40 


357.5 


34.8 


1.4 


0.0295 


30.250 


— > 


303.65 


2.8 


34.7 


1.1 


0.0259 


30.121 


18 


2005 TN53 


289.14 


84.8 


9.3 


25.0 


0.0652 


30.213 


— > 


250.71 


88.7 


9.3 


25.1 


0.0652 


30.099 


33 


2005 T074 


270.02 


301.7 


169.4 


5.2 


0.0501 


30.222 


— > 


230.28 


306.9 


169.4 


5.2 


0.0517 


30.095 


37 


2006 RJ103 


242.91 


24.2 


120.8 


8.2 


0.0276 


30.112 


— > 


202.77 


29.5 


120.9 


8.2 


0.0318 


29.990 


22 


2007 VL305 


354.05 


215.3 


188.6 


28.1 


0.0648 


30.081 


— > 


317.23 


217.5 


188.6 


28.2 


0.0620 


29.982 


20 



(Figs. 2 & 3) which are made for artificial Trojans with ini- 
tial eccentricity eo ~ 0. Since they are all around the leading 
Lagrange point (La), we should put them on Fig. 2 in Paper 
I (for Li) to compare. A comparison tell clearly that five 
of them are inside the most stable regions (Regions A and 
B) and the only exception is 2001 QR322 (hereafter QR322 
for short), whose position in the dynamical map is near the 
low-inclination end of the "arc structure" . We have shown 
in Paper I that this arc structure is related to a resonance 
as fazi — 2/ CT = ge, i.e., CI in Eq. (6) of this paper. 

A more accurate comparison can be performed by us- 
ing the dynamical maps at different inclinations in this paper 
(e.g. Fig. 3) . Recalling that all dynamical maps in this paper 
are for orbits around the trailing Lagrange point (L5), we 
cannot put the observed orbits directly on the correspond- 
ing maps. Nevertheless, since the absolute symmetry be- 
tween I/4 and Ls have been proven (Nesvorny & Dones 2002; 
Marzari, Tricarico & Scholl 2003a; Zhou, Dvorak & Sun 
2009a), it is possible to find the symmetrical counterpart 
around L5 point for each orbit. To do so, we read from Fig. 2 
in Paper I the libration amplitude Act for each observed 
NT, which have been listed in the last column of Table 1. 
Then, remembering the inclination and neglecting the semi- 
major axis, the corresponding counterpart around L5 is the 
one librating with the same Act value on the (00, eo) — plane 
with the same inclination. For example, the object QR322 
has a small inclination (1.3°), and approximately we regard 
it as on the slice of irj = 0°, i.e. the first panel in Fig. 3. 
We check the libration (Act = 77°) and find that the best 
counterpart of QR322 is (o ,eo) ~ (30.38 au, 0.290). It is 
on the stable segment to the right of the unstable gap at 
ao ~ 30.375 au. Similarly, we can locate the symmetrical 
correspondence of 2005 T074 (hereafter T074) in Fig. 5 for 
io = 5° at (ao, eo) ~ (30.26 au, 0.052). It is on the stable side 
of the boundary between stable and chaotic motion. More- 
over, it's in the very vicinity of the S3 resonance curve. We 
may conclude that the orbits of both QR322 and T074 are 
stable and these objects, however, are probably influenced 
by resonances like CI in Eq. (6) and S3 in Eq. (7), respec- 
tively. The rest four objects listed in Table 1, again this time, 
are well inside the most stable regions. Thus we neglect here 
the comparisons in detail. 

However, the observed orbits are not exactly on the 
slices we have calculated. To confirm our conclusion and to 
get more accurate estimation of the orbital stabilities, we in- 
tegrate one hundred clones of each nominal orbit (note the 
numbers listed in Table 1 have been rounded for simplic- 
ity) and check their stabilities using the frequency analysis 




a(AU) 

Figure 8. The orbital stabilities of 100 clones of 2006 RJ103. The 
color indicates the spectral number that is used as the stability 
indicator. The black stars designate the nominal orbital elements 
on the (a, e) and (a, i) plane (upper and lower panel respectively). 

method adopted in our papers. Considering the uncertain- 
ties from observation and orbital determination, the clone 
orbits are generated using the covariance matrix given in 
the AstDyS webpage. We thank Dr. Antonio Giorgilli for 
supplying us the necessary computing codes for generating 
the proper initial conditions. 

The results confirm the orbital stabilities of 2004 UP10, 
2005 TN53, 2006 RJ103 and 2007 VL305. They are deeply 
inside the stable region, as proven by the fact that nearly all 
the clones of them have SN< 50. As example, we present in 
Fig. 8 the SN of clones of 2006 RJ103. Out of the 100 orbits, 
only 4 have SN> 50 with the largest SN= 57. 

As mentioned above, T074 is close to the edge of stable 
region. Considerable fraction (~ 50%) of the clones has rel- 
atively large SN~ 70. And there is a trend of declining SN 
as the eccentricity of clones decreases. Thus we may expect 
that further observations will constrain the eccentricity to a 
smaller value. 

The most interesting result is from the case for QR322, 
as shown in Fig. 9. Because QR322 is in a narrow stable 
"stripe-like" area separated from the main stable region in 
the dynamical map, some clones of this object are protected 
by the C-type resonance /2a — 2/ CT = ge, while some other 
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Figure 9. The same as Fig. 8, but for object 2001 QR322. 



clones nearby may suffer the destroying effects from reso- 
nance like 2/ CT — /a:i — g + g& + gs — and its analogues. 
In Fig. 9, we see the stability of clones depends sensitively 
on the semimajor axis, which can be explained by the sensi- 
tive dependence of the libration frequency f a on semimajor 
axis. The eccentricity and inclination within the observing 
error range, on the other hand, do not influence the orbital 
stability considerably. Similar results can also be found in 
(Horner & Lykawka 2010b), where the lifetime of clones has 
a sharp edge on semimajor axis (Figs. 3&4 therein). In this 
sense, the QR322 is not a typical (thus not a good) example 
of NTs. If using it as a seed to generate a cloud of arti- 
ficial asteroids around to represent the real NTs (Horner 
& Lykawka 2010a), additional attention should be paid to 
make the conclusion drawn from the orbital simulations of 
the cloud consistent and convictive. 

In addition, we have integrated the six nominal orbits 
up to the Solar system age (4.5 Gyr) and found that all of 
them survive on the Trojan-like orbits. According to our 
investigations above, it is an expectable outcome for five 
objects. But for QR322, we would like to say, it is more or 
less by luck. Not like others, only a very slight deviation 
from the nominal orbit of QR322 may lead to an unstable 
orbit. 

In Fig. 10, we show such an example. The initial orbit 
is exactly the literal one listed in Table 1, in which the num- 
bers have been rounded for simplicity, e.g. the inclination 
1.3° is rounded from the nominal value 1.3220°. At the be- 
ginning 1 x 10 yrs, the orbit has a quite regular behaviour. 
The eccentricity is small, the inclination is constrained be- 
low 4 c irc, and the resonant angle a varies with an amplitude 
smaller than 70°. In most of the time, it is trapped in the 
i/i8 resonance, as AQ librates around 0° with an amplitude 
smaller than 180°. Then from 1.05 x 10 yrs, the eccentric- 
ity increases and the orbit loses stability, escaping from the 
Trojan-like orbit. 

A close observation at Fig. 10 may reveal that as early 
as 4 x 10 4 yrs, some chaotic features can be found in the 
semimajor axis' behaviour, and the eccentricity begins to 
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Figure 10. The orbital evolution of a clone of 2001 QR322. From 
top to bottom, we illustrate in 8 panels the temporal evolution 
of a,e,i,ui (critical angle of Kozai resonance), Q,a = A — Ag 
(resonant angle of the 1:1 mean motion resonance) , Ara = ra — rag 
(critical angle of vg resonance) and Af2 = f2 — f2g ( critical angle 
of fi8 resonance). 



show a slow increasing trend from 8 x 10 4 yrs. We think 
it is due to the fact that this orbit is in the resonance 
2/ CT — /2:i — g + <?g + ge — 0. The libration f a and apsi- 
dal precession g axe involved in this resonance, therefore the 
chaos was introduced to the behaviour of a and e. Finally, 
the resonance results in totally chaotic orbit. 



6 CONCLUSIONS 

As a continuation of our research in a previous paper on 
the dynamics of NTs, we present in this paper a detailed 
investigation on the stability of orbits on the initial plane of 
(ao,eo). By constructing dynamical maps on slices with dif- 
ferent inclinations io , we obtain a global view of the stability 
of NTs in the whole orbital parameter space. 

In the three most stable regions in inclination, we found 
in the representative dynamical maps the extension of sta- 
ble motion in eccentricity is 0.10 for low-inclination orbits 
with io < 12°, 0.12 for medium inclination 22° < io < 36°, 
and 0.04 for high inclination between 51° and 59°. The fea- 
ture of dynamical map was enriched by fine structures in it, 
indicating the diversity of orbital behaviour of NTs. 

Using the dynamical spectrum method based on the 
frequency analysis, we figured out the mechanisms creating 
the fine structures in the dynamical maps. We found two 
types of resonance may involve deeply in the dynamics of 
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NT. One is the secondary resonance (C-type) concerning 
mainly two frequencies: the libration frequency of the reso- 
nant angle a = A — Ag and the frequency of the quasi 2:1 
MMR between Neptune and Uranus 2Ag — A7. The other 
is the secular resonance in general sense (S-type), charac- 
terized by the commeasurability between different combina- 
tions of secular frequencies, such as the frequencies of apsi- 
dal and/or nodal precessions of NT and/or planets. Among 
the well-known secular resonances, the v$ resonance is very 
effective in driving a Trojan to highly eccentric orbit thus 
make a deep unstable gap in inclination around 44° . On the 
contrary, the z^is resonance, found nearly everywhere at low 
inclination, is so weak that it has hardly any influence on 
the dynamics of NT. 

Our study is not just theoretical, since we can place 
the observed NTs on the dynamical map and check whether 
they are trapped in or close to some identified resonances. In 
this way, the dynamical features of the observed objects can 
be predicted. We found 2004 UP10, 2005 TN53, 2006 RJ103 
and 2007 VL305 are deeply inside the region of the most 
regular motion, far away from dominant destructive secu- 
lar resonances. Therefore they must be stable. The orbits of 
2001 QR322 and 2005 T074 may survive to the Solar sys- 
tem age, but they are probably influenced by specific C-type 
and S-type resonances, respectively. And future observations 
may constrain their orbits to stable region a little deeper. 
To confirm the above conclusions, hundreds of orbital clones 
of the six observed objects were generated within the obser- 
vational error bars, their orbits were integrated numerically, 
and their stabilities were carefully inspected. 

Surely, the resonances listed in Eqs. (6), (7), (11), (12), 
(16) are not all of the possible resonances involved in the 
orbital evolution. Especially, we may miss some important 
resonances with higher order. However, even only these res- 
onances listed in our paper, have built a "resonant net" 
that connect different parts of the whole orbital parame- 
ter space. Correspondingly, there must be such a net in 
the phase space. With the help of this net, an orbit may 
diffuse in the phase space, resulting in significant change 
in the orbital behaviour. Long-term diffusion may be ex- 
pectable along the net. The most attractive diffusion, prob- 
ably is the slow transferring of an orbit from low inclina- 
tion to high inclination. Although in our preliminary tests, 
we did not find such inclination-increasing phenomenon, we 
would argue that this is a very important issue. The impor- 
tance arises not only from the interesting dynamics itself, 
but also from the fact that this problem is closely related to 
the origins of those highly inclined NTs. Again, their origins 
are closely related to the evolution of our planetary system 
in the early stage (Kortenkamp, Malhotra & Michtchenko 
2004; Nesvorny & Vokrouhlicky 2009). 

Another possible hint about the early evolution of the 
Solar system is buried in the C-type resonances. On one 
hand, C-type resonances are quite powerful in featuring the 
dynamics of NTs. On the other hand, these resonances are 
sensitively affected by the frequency /2:i. If the planetary or- 
bital configuration changed, even only very little, f-2-.i would 
change and this would lead to significant varying in the po- 
sition and strength of the C-type resonances. Thus a careful 
discussion on the influences of the C-type resonances, un- 
der current planetary configuration and under other possi- 
ble configuration, is critical. Similar analysis have been per- 



formed in the case of the formation of the Hecuba Gap in 
the main asteroid belt, where the Great Inequality (quasi 
5:2 MMR between Jupiter and Saturn) plays the role of /2a 
here (Henrard 1997). A closer analogue is about Jupiter Tro- 
jans (Robutel & Bodossian 2009) . For NTs, some work have 
been done (Nesvorny & Vokrouhlicky 2009; Lykawka et al 
2010), and our contribution is also undergoing. 
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